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Abstract 

We develop an iterative technique for computing the unstable and stable eigenfunctions of 
the invariant tori of difFeomorphisms. Using the approach of Jorba, the linearized equations are 
rewritten as a generalized eigenvalue problem. Casting the system in this light allows us to take 
advantage of the speed of eigenvalue solvers and create an efficient method for finding the first 
order approximations to the invariant manifolds of the torus. We present a numerical scheme 
based on the power method that can be used to determine the behavior normal to such tori, and 
give some examples of the application of the method. We confirm the qualitative conclusions of 
the Melnikov calculations of Lomeli and Meiss (2003) for a volume-preserving mapping. 

1 Introduction 

Hyperbolic sets commonly occur in dynamical systems; consequently, developing an 
understanding of the structure of their invariant manifolds is of considerable signifi- 
cance. For example, the stable and unstable invariant manifolds of a hyperbolic fixed 
point often intersect transversally, and the resulting heteroclinic tangle is a primary 
mechanism for the onset of chaos [GH83, Wig94 . The theory of transport in dynam- 
ical systems is based on the construction of regions bounded by partial barriers that 
can often be built using segments of stable and unstable manifolds [M ei92l |Wig92| . 
Moreover, Arnold's mechanism for drift in near integrable Hamiltonian systems is 
also based on the construction of transition chains made up of invariant manifolds 
|Arn64^ ILoc9 9J . In this paper we will develop a simple iterative method for the com- 
putation of such manifolds that was suggested to us by Hector Lomeli. 

The easiest algorithms for computing invariant manifolds of a hyperbolic set for a discrete 
dynamical system / are based on iteration. For example, suppose that / has a hyperbolic fixed 
point p. Near p, W'^{p) can be approximated by the linearized unstable space, E^{p), or more 
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accurately as the fixed point of a contraction map on sequences of points in a neighborhood of p 
|HDVS03] . Given such an approximation, one can start the iteration with a fundamental domain 
U on W^ip), i.e. a shell enclosing p whose inner boundary iterates to its outer boundary |LM00j . 
Iteration of U will then cover W^{p). The shell U can be covered with a grid of points on several 
spheres. Since distances between grid points typically grow exponentially with iteration, it is 
worthwhile to use distance as a diagnostic for inserting or removing additional grid points as 
needed. This method works well for one-dimensional manifolds Hob93j , and for two-dimensional 
manifolds when the two multipliers have the same magnitude _LM98. iLMOOj . More sophisticated 
methods for constructing invariant manifolds define a flow on the manifold that does not suffer 
from the problems of exponential nonuniformities when the multipliers have distinct magnitudes 
|K()98a,llK098hll(;vn4j . 

More generally, let / : M" — > M" be a diffeomorphism that defines the discrete dynamical system 

y' = f{y) , (1) 

and suppose that / has a compact invariant set T. Given a point z £ T, the linearization of the 
dynamics near T defines a skew-product system on (x, z) € M" x T 

xt+i = Df{zt)xt , 

zt+i = f{zt) ■ (2) 

This is an example of a cocycle over the dynamical system f\j- ^KH 99 , Supplement]. Many prop- 
erties of such systems are well-known. 

We will specialize to the case when T is an embedded invariant torus. Invariant tori commonly 
occur in dynamical systems. For example, KAM theory implies that nearly integrable, symplectic 
twist maps have Cantor sets of tori. Tori whose dimensions are less than the number of degrees of 
freedom of the system can be normally hyperbolic |Tre94l iBTOOj . Hyperbolic, or "whiskered" tori 
provide a primary mechanism for transport in these systems and are a necessary part of Arnold's 
construction of a transition chain for systems with more than two degrees of freedom. Invariant 
tori are also a common feature of volume-preserving mappings |PF88| IFKP89| IHZD98| IMLW98( 
ICMPTO2I ILMOOl ILM03j . Our goal in this paper is to develop an iterative algorithm for constructing 
the linear stable and unstable spaces of an invariant torus. 

Suppose that the function z : ^ T C M"" represents the embedding of the invariant r- 
dimensional torus T. We will use angle coordinates 9 on T*"; they are taken mod 1, so that 
z(6 + m) = z{6) for all m € Z'". We assume that the dynamics on T is conjugate to the rigid 
rotation 6 ^ 9 + iv with rotation vector G M*" so that 

fizi9)) = z{9 + u;). (3) 

For this case the cocycle ^ can be written 

x' = A{9)x , 

9' = 9 + L0, (4) 

where A : T'' ^ GL{n) is defined by A{9) = Df{z{9)). 
We will always assume that uj is incommensurate: 

co-m^n M meT'\ {0}, n e Z . (5) 
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Denoting the distance to the nearest integer lattice point as 



I \a\ \jk = min | |a — n| | , (6) 

for an a G M^'/ then a rotation vector oj is incommensurate if ||m • is nonzero whenever m is 
nonzero. When oj is incommensurate, the system Q) is a linear, quasiperiodic skew-product. 

Often it will also be necessary to assume that the rotation vector lo satisfies a Diophantine 
condition, i.e., that there exist c > and r > r such that 

||m-w||z >c/||m|r, VmGZ'^\{0} . (7) 

For example, KAM tori are Diophantine. As we will see next, many of the theorems that show @ 
can be simplified, or "reduced," rely on such conditions as well. We will often use the golden mean 

7 = ^—, (8) 
which is Diophantine and has the maximal value for c. 



2 Reducibility 

Following Lyapunov, a skew-product is reducible if there is a continuous change of variables C : 
T*" — > GL{n), so that defining x = C{9)y, @ becomes 

y' = Jy 

6' = e + uj , (9) 

such that the matrix 

J = {c{e + uj))-^A{9)c{e) (10) 

does not depend on 9 |Eli98llJor01j . The reducing transformation is the quasiperiodic generalization 
of the Floquet representation in the periodic case. Just like in the Floquet case, the transformation 
C and the matrix J are generally complex; however, they can be made real at the expense of 
allowing C to act on some finite cover of T'' |.TS8H IE.T82j . 

The main point of reducing the matrix is that the fundamental matrix solution of Q), 

t-i 

A('H9)^llA{9 + ju;), (11) 

j=0 

then has the representation 

A^'\9) = C{9 + tuj)j\C{9))-^ . 

Thus the stable, unstable, and center spaces of are easily found from those of the constant 
matrix J. Since the dynamics are now are trivial, the case that @ is reducible is ideal. However 
in contrast with the Floquet case, a reducing transformation cannot always be found. When A is 
smooth and close enough to a parameter-dependent, constant matrix and u) is Diophantine (O, then 

^Note that this expression is not a norm. 
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A is reducible for all but a small, positive-measured set of parameter values jJS92j . On the other 
hand, some systems with strongly varying A, like the skew product version of the quasiperiodic 
Schrodinger equation 

e{Un+l + Un-l) + V{6)Un = EUn 

with |e| < eQ{V,oj) are not necessarily reducible pioTl IFS WQOl IeIIQT] . 

One way to try to find a reducing transformation C{0) is to find solutions to the eigenvalue 
problem |Jor01j 

A{e)^{d) = XT^^{d), (12) 

where is the translation operator, \T^'il^]{0) = i{;{9 + uj). The goal is to find a set of eigenfunctions 
ip G C(T^, C") \ {0} and their corresponding eigenvalues A G C. It is easy to see that if {X,ip{6)) 
is an eigenvalue-eigenfunction pair for dJ, then for any m G Z^, 

is as well. A pair of eigenvalues such that A2 = Aie^'^*™ '^ are called "w-related." 

Whenever (jlj) is reducible then eigenvectors for (fT^ exist; indeed, it is not hard to prove the 
following result. 

Lemma 1 f | Jor01| ). Suppose A G C(T'^ ,GL{n)) is reducible to the constant matrix J through 
nU) . Then 

• if J has an eigenvalue-eigenvector pair {X,v) of J, the eigenvalue problem has a solution 
{X,C{e)v), and 

• if X is an eigenvalue of Ill^A) . then 3m £ II such that ^g^'^*'"''^ is an eigenvalue of J. 

Conversely, suppose there exist eigenfunctions ipj{6) G C(T'",C"'), j = l,...,n for il^) that are 
independent at each value of 9. Then A is reducible to the constant matrix J = diag{Xi, A2, • • • , A„) 
with C{6) = [V-i, V'2, • • • , V'n]- 

Note that w-related eigenfunctions are linearly dependent for each 9, even though as functions on 
the torus they are independent over C". Since the set {exp(— i27rA; ■ 9) : k € Z^'} is a Fourier basis 
for C(T'",C), when there is a set of eigenfunctions ipj{9) that are independent for each 9, then the 
set of functions 

{e-'^^^-'^j{9):j = l,...,n,k^r-} 

is a basis for C(r,C"). 

In addition to the right eigenfunctions, (|12j) . we will use left eigenfunctions, (j) G C(T'", C") \{0}, 
which are defined to be solutions of 

A*{9)T^<P{9) = Xm , (13) 

where A* = is the Hermitian conjugate of A. The left and right eigenfunctions (if they exist) 
are orthogonal for each 9. To prove this, we first state a useful lemma (see, e.g. |Los89) l. 

Lemma 2. Suppose g G C(T'",C), and 

g{9) = ag{9 + oj) (14) 

for incommensurate uj and some constant a G C. Then either g = or there is an m £ II 

and a constant c such that a = g"^'^*™''^ and g = ce^'^*™'^. 
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Proof. Expand g in a Fourier series on the torus with coefficients gm, m £ Z"^. If g ^ then there is 
least one m such that g^n 7^ 0. The Fourier transform of (|14() imphes that gm = ae^'^^^"'''^ gm-, so that 
a = e"-^'^*"*''^. Suppose there is some other nonzero Fourier coefficient, say gi, then ae^"^"^^'^ = 1, 
which implies {I — m) ■ oj £ "L^ but by © the only solution of this is / = m. Thus g has a single 
nonzero Fourier coefficient. □ 

Lemma 3. If ui is incommensurate, the right hl^) and left eigenf unctions corresponding to 
ijO -unrelated eigenvalues are orthogonal for each 0. 

Proof. Multiply <i)*{6+u:)A{e) = Xj(j)*{6) by ipkiO) on the right, and multiply A{9)^k{0) = >^ki^k{6 + 
uj) by + ^) on the left. Then subtract the results to obtain 

If j = k, then this equation implies that (l)*j^{9)'iljk{0) is constant on a dense set (since uj is incom- 
mensurate) and so must be constant everywhere by continuity. The constant can be chosen to be 
one without loss of generality. If j 7^ k, then Lem. [51 implies either that Xj = A^e"^'^*™''^, which 
would mean that the two eigenvalues are cj-related, or (j)*{6)ipk{G) = 0. □ 

Just as a constant matrix need not be diagonalizable, a complete set of eigenfunctions of A need 
not exist, even when A is reducible. To deal with this case, define a generalized eigenfunction for 
eigenvalue A to be a nonzero, continuous solution to 

\{[A{e + tu:)-\T^I\^{d)=^ , (15) 
t=o 

for some A; € N. Suppose that for a fixed A, there exist exactly / continuous solutions ipj,j = 1, . . . , Z 
to (|15j) that are independent for each 9. These are a basis for an /-dimensional vector space 
Ex C C(T''",C"). It is easy to see that E\ is an invariant subspace under A. Indeed, if -f/; G E\, 
then 

{A{e-uj)-\T^i)^{e-oj) 

is in E\ as well since it is annihilated by H15|) . 

Just as there is a correspondence between regular eigenfunctions of A and eigenvectors of J, 
there is also a correspondence between generalized eigenfunctions and the generalized eigenvectors 
of J. 

Lemma 4. The matrix A{9) is reducible if and only if there exists a set of generalized eigenfunctions 
ipj G C(T^', C"), j = 1, . . . ,n that are linearly independent for each 6. 

Proof. Suppose that A is reducible to J by the transformation C{9). If li; is a generalized eigenvector 
of J with eigenvalue A, then for some G N, ( J - Xir w = 0. Defining ^{9) = C{9)w, 

(J -XI)w = [C{9 + uj))-^C{9 + cj) (J - XI) C-\9)Ci9)w 
= {C{9 + uj))-' [A{9) - XT^I]i^{9) . 

Applying (J— A/) an additional k—1 times shows that ip satisfies (|15|) . Since C{9) is invertible, linear 
independence of the generalized eigenvectors implies the linear independence of the corresponding 
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To prove the converse, let E = span(V'i, V'2, • • • , fpn) C C(T^, C"). Since is a finite-dimensional 
vector space, the standard representation result for generalized eigenfunctions |HS74[ Appendix 
III] implies that E is the direct sum of generalized eigenspaces Ex., j = corresponding to 



eigenvalues Xj, and that each of these has a basis ipj ,k = 1, . . .nj = dimE\. such that 



[A{9) - XjT^] iP^+\9) =7p^{e + Lo) , k = l...nj-l. 

Define the matrix C{6) = [ipl{6), . . . (0)]. By assumption, detC{6) 7^ 0. Then by the reverse 
of the previous argument, it is easy to see that J{6) = (C{9 + uj))^^A(6)C{6) does not depend on 
9, and has generalized eigenvectors Wi = {C{9))~^il)i{9). □ 

In this paper, the skew-product system (jlj is the linearization of a map on an embedded 
invariant torus. For this case there are always unit eigenvalues: 

Lemma 5. /// has an embedded r -dimensional invariant torus on which the motion is conjugate 
to a rigid rotation, then the eigenvalue problem US]) associated with the linearization ^ has r 
linearly-independent eigenfunctions ^i{9) = ^(9), i = 1, . . . ,r with A = 1. 

Proof. Differentiation of the conjugacy (jSJ with respect to 9i shows that the functions ^/Ji are 
eigenfunctions with eigenvalue one. Since z{9) is an embedding, Dqz has rank r for all 9 in T'', so 
the eigenfunctions are independent. □ 

The reducibility of A thus hinges on the existence of n — r generalized eigenfunctions that are 
transverse to the tangent space of the invariant torus. When u; is incommensurate, each eigenvalue 
A is associated with a dense set of w-related eigenvalues. Thus the closure of the spectrum of 1)12^ 
will correspond to the union of up to n — r + 1 circles centered at the origin of the complex plane. 
As an illustration. Fig. ^shows the an approximate spectrum for a three-dimensional mapping (this 
mapping is studied in ^ with an invariant circle. For this figure, the eigenvalues were computed 
using a Fourier series algorithm developed by Jorba |Jor01j : only finitely many eigenvalues are 
found due to truncation of the Fourier series. 

The eigenvalues of (fT^ are related to the Lyapunov multipliers of the skew- product Q. Ac- 
cording to Oseledec's theorem, can be decomposed into invariant subspaces spaces Hi{9), i = 
1 . . . k{9) for almost all 9 such that for each xq G Hi{9), the Lyapunov multiplier, 

//i(xo) = lim ||^W(0)xo||^/* , 

exists. Moreover the matrix 

l/2t 



A{9) = lim \(a^^\9)X A^*\9) 



exists for almost all 9. An eigenvector v{9) of A is called a Lyapunov vector |.TOn4j . The corre- 
sponding Lyapunov multiplier /i is the growth rate of v under iteration 

^Ji{v) = lim \\A^''\9)v{9)f . (16) 

t— »oo 

In the "nicest" case, T is uniformly hyperbolic, and then the spectral subspaces Hi{z) vary 
continuously and the corresponding the Lyapunov multipliers are constant and exist everywhere. 
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Figure 1: Numerically generated spectrum of for the volume- preserving example 1281 1 using 16 complex 
Fourier modes. Here the parameters are e = 0, v = .35, and a; = 7" , the inverse of © . There are 99 = 3(2 x 16 + 1) 
eigenvalues shown; all but four fall very nearly one of three circles whose radii correspond to the three eigenvalues 
A = . The final four values are are numerical artifacts.. eps 



When A is reducible, the hmits exist uniformly for all 9 so that the matrix A is continuous. In 
this case the Lyapunov spectrum coincides with the magnitudes of the eigenvalues Aj. A slightly 
weaker converse is also true: 

Theorem 6 ([Eli98, JS8IJ). If A{6) has a uniform and simple Lyapunov spectrum > ^2 > 
. . . /i„ and if to satisfies a Diophantine condition then A is reducible. 

We will use this result in the next section to compute the eigenfunctions. 

The results in this section reveal the connections between the reducibility of a skew product and 
generalized eigenvectors: the transformation that reduces the linear quasiperiodic skew-product Q 
is constructed from n linearly independent generalized eigenfunctions (|15() and their corresponding 
w-unrelated eigenvalues. These functions will form the transformation matrix C{6). 

Our goal is to use the eigenfunctions as the linear approximations to the invariant manifolds 
Wiz^e)), W{z{e)), and W'{z{e)). 

3 Computing Eigenfunctions 

One method for finding the eigenfunctions of (|12j) is to expand T^.,A{9) and il^{9) into Fourier 
series and use the fact that the Fourier coefficients of a product are convolutions of the coefficients 
of the individual functions. If the series are truncated, then the resulting equations correspond 
to a finite dimensional linear system, and they can be solved by standard numerical techniques 
|Jor01j . Iterative refinement methods have also been developed that use Fourier series ^HDL L04] . 
This method has several drawbacks. First of all, the construction of the representation of (|12jl in 
a real Fourier basis is somewhat involved. Moreover, the dimensionality of the resulting operators 
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A(6) and can be extremely large, especially if r is large. A representation with N complex 
Fourier modes and their conjugates on T** requires (2A^ + l)*" coefficients. Since each component of 
ip requires such a representation, the truncated representation will be of length n{2N + 1)^. In our 
experience, a minimum of = 32 modes is required for reasonable (twelve digit) accuracy even in 
very smooth examples.^ Thus, if n = 3, and the invariant set is a circle, the representation of A 
requires a matrix of size 195 x 195. If the space were five-dimensional and the invariant set where 
a three torus, the matrices would have a minimum size of the order of 10^ x 10^. Even the most 
reliable eigenvector solvers could prove unwieldy in such a situation. 

Thus we are motivated to develop an iterative method for finding eigenfunctions. 



3.1 Power Method 



When the eigenvalues of a constant matrix have distinct magnitudes, the power method can be an 
effective technique for finding the eigenvectors. We will use this idea here to find the eigenfunctions 
of A{9) transverse to the torus (recall that the remaining r tangent eigenfunctions are known by 
Lem. 121). While the power method will find the dominant eigenfunction as long as A is reducible 
and the dominant Lyapunov multiplier is distinct, for the remainder of this section we will assume 
that A{9) is a real matrix with a uniform and simple transverse Lyapunov spectrum. Thus the 
transverse Lyapunov multipliers have distinct magnitudes: 



Ml > > > fJ-n- 



1 



. n 



If we also assume that uj satisfies a Diophantine condition , then Thm. El applies and A is re- 
ducible. Since the multipliers are distinct and A is assumed to be real, there will be a corresponding 
set of real eigenvalues Aj = ±/Xj and real eigenfunctions ipi of (|12j) . Algorithms for the more general 
case of multiple and complex eigenvalues will be discussed in a future paper |WM05j . 

As in the power method, we begin with an arbitrary initial vector q^'^\ Since A is assumed 
reducible, every vector can written as a linear combination of the eigenfunctions at some 



,(0) 



(17) 



and generically ai ^ 0. Iteratively define sequences u^''^ and q^^^ by 



Jk) 



and = ^(^0 + kuj)u^''^ , 



(18) 



This implies Hj^o H^^^^H = ||^^''^(6'o)9^°^||, where A^^) is defined in Since ipi is the eigenfunc- 
tion with the largest eigenvalue, we have 



u 



(k) 



||^W(0o)9(°)| 



(19) 



Because u^^^ is a unit vector by definition and ■i/'i is continuous, its coefficient in H19|) must be 
bounded. Defining Sk = sgn(A)'^, we therefore have 



^For example, /e in fTH with e < 0.35. 



L(^o+Mir 
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which is the vector we were seeking (up to the choice in sign). Moreover, the Lyapunov multipher 
is obtained simply from 

l/k 



k 




= lim 

fe— >oo 

however, this process has shortcomings in certain appUcations that the following modifications help 
to alleviate. 

To find "01(^0) we select a subsequence kj on which ||/cja;||zr approaches zero (using the pseudo- 
norm ©). For the case of an invariant circle, r = 1, the sequence kj can be constructed from the 
continued fraction expansion 

w = ao + 1/ (ai + 1/(02 + ...)) = [ao, ai, 02, • • •] , 

with elements Oj G for example, 7"^ = [0, 1, 1, 1, . . .]. An optimal sequence ki corresponds to 
the denominators of the continued fraction convergents, and are given by the recursion relation 

fcj+i = ajkj + kj-i , j e Z+ , ko = 1 ,ki = I . (20) 

For 7~^, this gives the Fibonacci sequence {1, 1, 2, 3, 5, . . .}. For higher dimensional cases, we could 
use the generalized Farey tree of |K086j . which also generates the best approximants. 

Given the sequence kj, a polynomial interpolation on u^^^^ through the nodes 6q + kjuj can 
be used to approximate the value of ipi at 60 . We stop the iteration once a desired accuracy is 
obtained. The Lyapunov multiplier fii can be rapidly and accurately computed by using V'i(^o) in 

m- 

If the system has a negative dominant eigenvalue, then the direction of the iterates will oscillate, 
and for this reason, we normalize the n^'^^^ so that one component is +1, say. If interpolation is 
carried out in this way, and an approximation to Tpi{6o) is accepted but future iterates of ^1(^0) 
oscillate in direction, then the eigenvalue must be negative. 

An alternative method for computing ipi{9o) is to begin with an arbitrary vector q^'^^ at a point 
9q — kuj and iterate it forward k steps to 6q. However, this does not give an efficient method for 
estimating the error, and if nothing is known about the system a priori, then it is difficult to select 
the appropriate k. Since, by (fT!l|) . the convergence rate of the iterates to -01 is A2/A1, the error 
in the interpolant sequence n^'^^^ should decrease geometrically whenever the Lyapunov spectrum 
is simple, and we use this as a convergence check. Finally, if the interpolation scheme does not 
converge, then this can be a indication that the system is either irreducible or has a dominant 
Lyapunov multiplier with multiplicity larger than one. Naively accepting a given iterate u^^'^ as the 
approximation to ^'(^o + ku^) would then be a grave mistake. 

Once we have obtained V'i(^o) and Ai to within a selected accuracy, we generate V'i(^) on a 
dense set of 6 values by defining 

+ ku:) = Ar*^A(')(0o)0i(^o) (21) 
The resulting function satisfies (|12j) since 

A(0o + fcw)V(0o + kuj) = A-^'^(^-+i)(0o)V'i(^o) 

= Xi;i9o + ik + l)uj). 
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When A is reducible the function ipi computed in this way should be continuous on the invariant 
torus. Under the assumption that the Lyapunov spectrum is uniform, the associated eigenspaces 
are themselves continuous. This implies the continuity of the direction of ipi{0), but its magnitude 
will only be continuous when the function c{9) defined by 

c(eo + M =Ar'P^')(^o)V'i(^o)|| (22) 

is continuous. While c{6) is bounded because HA*^'^) (0)^(0) || ~ A*"', its continuity is not completely 
obvious. We will see, however, that c will be continuous in the examples presented in 311 

3.2 Subdominant Eigenfunctions 

The power method easily finds the dominant right eigenfunction and its corresponding eigenvalue. 
It is also easy to find the eigenfunction with the smallest multiplier by iterating with {A{9))~^ and 
replacing u with —uj. Similarly, the dominant left eigenfunction can be found using the power 
method for (|13j) . or equivalently by using (|12|) upon replacing A{6) with A* (6), and iv with —uj. 
However, to find intermediate eigenfunctions, another strategy must be used. Here we discuss two 
possibilities. 

One such technique is a generalization of Hotelling's deflation method Ral65 . In this method 
we modify A{9) to create a new cocycle in which Ai has been replaced with zero. The power 
method using this modified skew-product will then converge to (A2,'02(^))! the next largest eigen- 
value/eigenfunction pair. The modified system can be constructed if we already know the dominant 
right and left eigenfunctions: 

Lemma 7. If A(9) is reducible and has dominant eigenvalue Ai with corresponding right and left 
eigenfunctions ipi{9) and </>i(0), normalized such that (/)l{9)ipi{9) = 1, then the matrix 

Wi{9) = A{9) - XiM^ + (23) 

has the same eigenvalues and eigenfunctions as A{9) except Ai has been replaced with 0. 

Proof. Since the left and right eigenfunctions corresponding to distinct eigenvalues are orthogonal 
by Lem. 131 the eigenfunctions of 1^1(6*) are the same as those of ^4(6*). To see Wi{9) and A{9) share 
the same spectrum except for Ai, apply {C{9 + <^))^^ on the left and C{9) on the right of (|23|) 
arriving at 

{C{9 + u;))-'Wi{9)C{9) = J - Xi[ei,0, - ■ ■ ,0] . 

□ 

Consequently, the power method can be used with Wi{9) to find ip2{9). Similarly, the corre- 
sponding left eigenfunction can be found using W^. This method can be continued in principle to 
compute all of the eigenfunctions using successive deflations 

k 

Wki9) = Aie) - J2 ^^MO + oj)(t>*{9), k = 2,...,n-l, 

1=1 

thereby determining 7p^{9), . . . , V'n(^)- 

A second method to find subdominant eigenfunctions is iterative Gram-Schmidt orthogonaliza- 
tion. In lieu of making a rank-one modification of A{9), this method removes any components in 
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the direction of the dominant eigenfunction at each iterate of the power method. Again, suppose 
that we have found the dominant right and left eigenfunctions. To determine 1^2, begin as usual 
with a randomly chosen q^^^ at an initial point 6. Since we assume A is reducible, q^^^ can be 
written in the eigenfunction basis as ()17() . Since ip2 must be orthogonal to we project out the 
■01 component by applying the transformation 

q--<l^ = q- JI[^^'M 9). (24) 

In principle, q± = ^"=2 '^jV'j(^) bas no component in the ipi direction and thus the power method 
applied to q± will converge to ip2- However, roundoff errors inevitably introduce a small component 
in the ■01 direction and so we reapply the orthogonalization after each iteration. Similarly 02 can 
be found by projecting out the 0i component and iterating with the adjoint matrix. 

The orthogonalization technique to find the remaining subdominant eigenfunctions requires 
projecting out the previously- found dominant components. 



3.3 Stability Analysis 

We first look at how errors in the computation of the dominant eigenfunctions affect the deflation 
process. Suppose that the dominant eigenfunctions have been obtained up to some error 

where C G M. Here we keep only the terms in the subspaces spanned by the two most dominant 
eigenfunctions and also assume the errors to be of the same order of magnitude. Then the deflated 
matrix Wi can be computed only approximately as 

Wi{0) = A{e) - XiMe + u)^i{9) . 

Now suppose qo = aiipi{9) + a2'02(^) is the initial guess in the computation of ip2- The error 
at the flrst step is 

(Wiie) - Wi{e))qo = Xie{aiMO + u;) + Ca2MS + ^)) + C(e^) = ^(Aie) 

Continuing this analysis shows that the error increases by a factor of roughly Ai for each iteration. 
Since the result at each step will be scaled by approximately A2 , the error in the second eigenfunction 
due to the error in ipi at the jth iterate will be 0(6 (A1/A2)''). This instability is also inherent in 
Hotelling's deflation method for constant matrices |Ral65j . 

The deflation method can be modifled to ameliorate this instability by applying the the projec- 
tion (|24|) to (approximately) remove the growing components in the tpi direction. This process can 
be applied occasionally as a corrective procedure. For eigenfunctions corresponding to eigenvalues 
of decreasing magnitude, this corrective procedure must be applied increasingly often in order to 
maintain the desired accuracy. Of course, the correction (|24|) can only be computed using approxi- 
mations to the eigenfunctions. Keeping this in mind, a straightforward calculation shows the error 
in the (p+ 1)** eigenfunction is of the same magnitude as the error made in the original p dominant 
eigenfunctions. 
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Note that Gram-Schmidt orthogonahzation amounts to applying ()24() at every iterate; therefore, 
Gram-Schmidt has similar stability properties to Hotelling's deflation with the projective refine- 
ment. So while the error in successive eigenfunctions for both methods is theoretically of the same 
order as the initial error in the dominant eigenfunctions, in practice the instabilities inherent in 
these two methods makes reliable computation of the subdominant eigenfunctions very difficult to 
achieve. Indeed, in |Ral65j it is remarked that the analogous deflationary procedures for constant 
matrices can often only be employed effectively by experienced numerical analysts. The deficiencies 
of the defiationary methods have prompted us to seek more stable techniques that will appear in a 
forthcoming paper |WMn5j . 



4 Examples 

As a first test of the power method, we construct a diffeomorphism on X T*" with a hyperbolic 
invariant torus whose manifolds are known explicitly. Using the coordinates {x, y, ^) G x T^', the 
map is defined by 

/l-6 + 5x + ey {Xg{e + uj) - 5g{e))\ 
fix,y,e)= Xy , (25) 

V d + LO I 

where g{0) = g{9 + m) \/m G 1/ . The torus T = {(1, 0, 0) : G T'^} is invariant with rotation vector 



ijj. 



Moreover, it is easy to see that the surfaces 



5 ^ {y = 0} 

U = {x = l + eyg{e)] 

are invariant under (|25|) : ifA>l>(5>0 these are the stable and unstable manifolds, respectively, 
of T . The mapping is volume and orientation preserving when A = 5^^. 
The linearization of H25|) on the torus gives the matrix 

(5 e{\g{9 + uj) - 5g{e)) 0^ 
A{e) = A 

VO /y 

For the case r = 1, T is an invariant circle, and the three w-unrelated eigenfunctions of ()12() . are 

V'l = 1 , ^^2 = , = , (26) 




with the eigenvalues A, 1, and 6, respectively. 

The power method of ^l-i. II works well to compute the eigenfunctions (|26|) . The convergence to 
V'l is governed by H19|l . which shows that the k^^ iterate of the approximate eigenfunction, u'^^\ is 
(up to normalization) 

u^^'^ oc aiV-i + X'''a2iJ2 + (jj "sV's- (27) 
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Equations ()26() and (|27() then imply that the error in the first component should decay like 1/A 
while the error in the third component should decrease as 5/\. Figure [21 is a log-linear plot of the 
error in the first and third components of the k^^ approximation to the dominant eigenfunction. 
As can be seen, the convergence rate is exactly as predicted. Redoing the calculation with various 
values 5 and A gives similar agreement with the theoretical prediction. 




20 k 40 60 



Figure 2: Error in the first (upper curve) and third (lower curve) components of the eigenfunction as a function 
of k, the number of iterates, for the test system 12511 . Here g{9) — 0.7cos^ + O.lSsinS, A — 1/5 = 2, e — .1, and 
LO = 7~^. Also shown are the predicted slopes of the errors. .eps 

To compute V'i(^o)> we interpolate as discussed in ^3.11 using continued fractions to generate the 
sequence (|2()|) of iterates kj. Generally, the error in a p^'^-order polynomial interpolation through 
nodes xq, . . . ,Xp for a C'P~^^ function at a point x is 

°(^n(— ))■ 

We chose a third-order interpolation scheme. Since the sequence of best approximants to uj obeys 
||%w||z = 0{l/kj), and the Fibonacci numbers grow as 7-', the third-order interpolation error 
should scale as 

This is confirmed by our computations, as shown in Fig. |2l Note that the asymptotic slope is 
achieved only after the transient terms lying in the subdominant eigendirections have been elimi- 
nated (the second and third terms in H27|) ). Fig. |31 demonstrates that more iterations are needed 
to remove this transient when the eigenvalue ratio A1/A2 = A is smaller. 

We accepted the value of ipi{0) when successive interpolations varied by less than 10~^^. Sub- 
sequently, we set 9 = 0, and v{0) = ipi{0) in (fT6|) in order to compute Ai. We accept the eigenvalue 
when the computations for the fcj and /cj+i iterations differ by less than 10~^^. The computed ipi{0) 
and Ai were then used in (|2H) to generate ipi^kuj) at 100 locations round the invariant circle. As 
shown in Fig. I^the sup- norm of the final error in tpi is less than 5 x 10~^^. The sup-norm of the 
error in the eigenfunction -03 — computed using iteration with — is of the order 10~^^. 

To test the method in a higher dimensional case, we set r = 2 in H25() . so that the invariant torus 
is two-dimensional. The rotation vector was taken to he uj = (r, ) , where r is the real solution 
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Figure 3: Interpolation error in the computation of tpi{0) for 12511 . For these curves, X = S ^, tj = 7, and for curve 
(A) A = 2, e = 0.1; (B) A = 1.5, e = 0.3; and (C) A = 1.1, e = 0.1. The predicted slope of -4 is also shown for 
reference., eps 




e 

Figure 4: Supnorm of the error in the numerical eigenvector tpi{9) for 12511 with the same parameters as Fig.|5|.eps 
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of the algebraic equation — x — 1 = 0. This rotation vector is Diophantine since the numbers 
(1, r, r^) form an integral basis for a cubic field |Cus72j . This time, instead of interpolating in order 
to find ■01 at a particular value of 6, we iterated 10^ times in order to remove any transient terms 
and simply accepted the last iterate as the approximation to ipi. We recorded the location 6*0 at 
which we have this approximation to tpi. We then set v{6q) = ipi{9o) in (|16|) in order to compute 
Ai, the value of which was accepted when successive approximations varied by less than 10~^^. 
Lastly, we generated iPi{Oq + ku) at 2 x 10^ locations on the torus. Fig. [HI is a plot of the first 
component of the computed -01- Because the locations Oq + ku) are irregularly spaced, Delaunay 
triangularization was used to compute the vertices of each triangular panel in the surface shown. 




Figure 5: Plot of the first component of ?/'i(^i,^2) for tlie map H25|l witli cu = (r, r^). The choice of parameters 
was e = .5, A = 1/5 = 1.5, and g{9) = cos(27r6li) cos(47r^2). The error in the computed values of ?/>! was less than 
5 X 10~^^.eps 

A more challenging test case is given by the volume-preserving example that appears in jLMOOl 
ILMD.Sj . Here, the unperturbed diffeomorphism /o on has a set of hyperbolic invariant circles 
whose invariant manifolds are also analytically known, but in this case they form heteroclinic 
connections between neighboring invariant circles. It was shown in |TMn3j using a Melnikov method 
that this connection is destroyed upon perturbation and the topology of the resulting intersections 
of the invariant manifolds undergoes bifurcations as the parameters vary. Our goal here is to explore 
the resulting manifolds numerically to verify and extend the predictions of the Melnikov theory. 

The unperturbed map /o in cylindrical coordinates (r, 6, z) is given by 

/ h-^{r + h{z)) -z \ 
h{r,e,z) = e + LO , (28) 

\ r + h{z) - 1 / 

where the function h{z) is an increasing circle diffeomorphism, i.e. h{z + 1) = h{z) + 1 Vz G M. For 
any such diffeomorphism, /q preserves the volume form dr Ad6 Adz. The map /o has an invariant 
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circle at T = {(1, 9, z*),9 G §} where z* is a fixed point of h. Tlie circle T has invariant manifolds 

U = {r = h-^(z) - h(z) + 1} , 
S = {r = l} . 

If A = TTj-^ > 1, then U is the unstable manifold of T and S is its stable manifold. Moreover, 

n'{z*) ' ' 

if h{z) has several fixed points, z* , then U and S intersect at each of the corresponding invariant 
circles at r = 1; thus, adjacent invariant circles have heteroclinic connections. It is interesting that 
with a general choice of diffeomorphism h, the map /o appears to be nonintegrable even though it 
exhibits heteroclinic connections |LM03j . 

Since the cylindrical axis r = and the set r > are invariant under /o, we can introduce 
rectangular coordinates {x,y,z) using the volume preserving transformation 

{r,e,z) {x,y,z) , 
X = V2r cos{2it9) , y =\/2r sin(27r0) , z = z , 

so that r = ^(x^ + y^) is the "symplectic" radius. We use the same notation for the map in the 
new coordinates, letting {x' , y', z') = fo{x, y, z). Now /o preserves the volume form dx A dy A dz. 

The linearization of fo{x,y,z) on the circle T = {x^ + y'^ = 2, z = z*} gives rise to a skew 
product of the form @ with the matrix 

^cos(27ra;) + 6 cos{2tt9) cos{2tt{9 + lu)) - sin(27ru;) + 6 sin(27r6') cos(27r(6' + lo)) 
Ao{9)= ( sin(27ru;) +(5cos(27r6l)sin(27r(6' + u;)) cos(27ra;) + 5sin(27r6l) sin(27r(6' + w)) 

^/2cos(27r6l) \/2 sin(27r6l) A^^^ 



(29) 



where 5 = X — 1. 

It is straightforward to verify that eigenvectors of H12|) for the matrix ()29|1 are 



/cos(27r6l)\ /- sin(27r6')\ 

V-i = sin(27r0) ,ij2=\ cos{27t9) , 

with eigenvalues Ai = A, A2 = 1, and A3 = A^^. It is also easy to verify that ipi is tangent to U, il)2 
to the circle T, and V's to 5. 

An especially simple case occurs when h is chosen to be the function |LMfl3j 

... 1 ^ / (1/ + 1) tanfvrz) + 1/ - 1 
hu{z) = — arctan 



TT \{v — 1) tan(7r2;) + v + 1^ 

In addition to being rotationally symmetric, /o now preserves the function 

J(r, 9, z) = 2v cos(27rr) + (1 - v^) cos(27rz) sin(27rr) ; (30) 

consequently with this choice of /i, /o is integrable. In this case the surfaces lA and S lie on the 
level set J = 2v. This particular choice for the circle diffeomorphism is especially nice because 
iteration is easy: /i^(-z) = hi^k{z) , VA; G Z. The function h,y has fixed points at ±z* = ib|, with 
h'(±z*) = z^^^. When u < 1, the set S is the stable manifold of the circle 71 at —z* with eigenvalue 
u, and the unstable manifold of the circle 7+ at +z* with eigenvalue 
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To destroy the heteroclinic connection between the circles 7±, we perturb /o by composing 
it with a near-identity transformation, Id + eP. The resulting map is still volume preserving 
for all e whenever the Jacobian DP is nilpotent, i.e., {DPi)^ = |LM00j . We select two such 
transformations, so that 

A = (/d + ePi)o(/d + eP2)o/o , (31) 

where 



P2ix,y,z) = (0,0,r-l) . 

Each of these Pi clearly have nilpotent Jacobians, and they also vanish on the invariant circles 
T±, so that no additional numerical work to find the circles is needed. Moreover, since the maps 
Id + ePi are diffeomorphisms, /e is a volume-preserving diffeomoprhism for all e. In this case the 
perturbed skew product has a matrix defined by 

/I - 2z*e'x'{l + {y'f) -2z*e^y'{l + {y' f) -2z*e{l + {y'f)\ 
A,{e) = 1 ^o(^) 

V ex' ey' 1 / 

where (x', y') = {^/2cos{2^T {9 + uj)), -v/2 sin(27r(0-|-u;)) corresponds to the rotated point on the circle. 

To compute an eigenfunction, we begin with an arbitrary initial guess go- Iteration is carried 
out using (|18)) . normalizing the iterates to have sup-norm one at each iterate. As before we use 
the sequence of convergents ki given by ()2U() and a third order polynomial interpolation to to find 
'i/'i(0). We continue iterating until the difference between successive interpolations differs by less 
than 10-12 

Once we find V'i(O)! we set f(0) = V'i(O) i^i (|lfc>|) and compute the corresponding eigenvalue. 
Convergence to the eigenvalue is declared when computations at the ki^ and /cj+i*'^ iterates vary 
by less than 10"^. Once A is known, we then produce V'i(^) s-t different locations on the torus using 
(|2H) . This creates a smooth eigenfunction in T, as can be seen in Fig. El here we show HV'iC^)!! 
and zoom-in near 9 ~ 0.955. A three-dimensional plot of ipi constructed from 400 iterates of ^(0) 
for the same parameters is shown in Fig. EJ Similarly, backward iteration with the power method 
can be used to find the stable eigenfunction. 

When e = the eigenvalues are (i^-^, 1, y). We show in Fig. IHlhow they vary with e for three 
different choices of rotation number. Note that the eigenvalues appear to approach one at critical 
values of e. Since the perturbation is chosen to preserve the circle and its rotation number, there is 
the slight possibility that these locations correspond to transcritical torus bifurcations. However, 
we have not been able to verify this as most trajectories near the tori are unbounded. Note that 
the power method becomes increasingly difficult to apply when the eigenvalues approach 1 because 
the linear convergence rate of the iterates u^^^ is given by A2/A1 ~ 1. While this fact may explain 
some of the gaps in the graphs in Fig. |H1 we are currently using more sophistocated techniques to 
determine exactly what happens dynamically in these regions. 

Once we have obtained the eigenfunctions ipi{9), we can plot the manifolds iteratively. Here we 
will compute the unstable manifold of 7+ and the stable manifold of 7^ . To draw the manifolds 
we first construct an approximate fundamental domain |LM00j . A fundamental domain for a two- 
dimensional invariant manifold is an annulus S with the property that the components, Cj, of its 
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Figure 6: The sup-norm of ip{9) for 71 generated using ll^ for the map with f = 0.55, uj = -y and e = 0.4. 
To show that the generated eigenfunction is continuous, we zoom-in on a selected point, .eps 




Figure 7: Plot of the unstable eigenfunction ipi on the circle T_ for the map 1311 1 with the same parameters as 
Fig.|^ .eps 
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Figure 8: Eigenvalues, Ai, of the circle T- for the map l|31|l as a function of the perturbation strength e with v — 0.55, 
and llI — 7"^ (red curve), 7""^ (green curve) and 7"* (blue curve) as a function of e. For e = Ai = X^^ — v. Since 
the map is volume preserving the product of its eigenvalues is one. .eps 
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boundary, dS = ci U C2, are related by iteration, /(ci) = C2- Using this fact, a fundamental domain 
can be iterated forward or backward to form an entire branch of the two-dimensional unstable or 
stable manifold of an invariant set. 

For the case of an unstable manifold, we form an approximate fundamental domain from ipi 
by choosing a factor /) <C 1 and defining a mesh of points from pX^^"^iJj{9j), k = 0, . . . ,m, j = 
0, . . . ,N so that the inner loop, ci, corresponds to the points with k = and the outer to A; = m. 
Quadrilateral panels are drawn using OpenGL to connect the points on this mesh, and then each 
point on this mesh is iterated forward to form a new annulus f{S) that is concatenated onto S. 
Panels are again drawn connecting the new grid of points on f{S) and the process is repeated. In 
our computations we choose p = 10^^, m = 5, and N = 400. 

In the mapping investigated here, the manifolds have transversal intersections as was predicted 
by the Melnikov calculations in |LM03j . For this reason, we cannot iterate a fundamental domain 
too many times before the oscillations of the manifolds become wild and the pictures too difficult 
to interpret. The number of iterates will be of order — log p/ log A for the manifolds to grow to a 
size of order 1. We used 14 iterates in the plots shown in Fig. ini^a) and (b), 12 for (c), and 27 for 
(d). 

In Fig. 1^ we show the invariant circles T± (green curves) as well as the unstable manifold of 7+ 
(red surface) and the stable manifold of 7L (blue surface). The first two panels, (a) and (b), show 
the manifolds intersecting along an infinite spiral that is forward asymptotic to 71 and backward 
asymptotic to 7^. In the last two panels, (c) and (d) the intersection curves undergo a reconnection. 

Since the images of each fundamental annulus S cover the invariant manifold, if there are inter- 
section curves, they must go through each fundamental annulus on each manifold. A fundamental 
annulus can be considered topologically to be a torus if we identify the boundary ci with C2 using 
the fact that they are images under /. Thus we can classify the intersection curves as elements of 
the homology group, 7?, of the torus. The intersection curves in Fig. EJa) and (b) have homology 
(1,0), crossing the fundamental domain once vertically with no net rotations in 9. As u; is increased 
to values near the resonance 1/3 there is a homology bifurcation. This is shown in Fig. |^ where 
Lo = ~ 0.382. Though the intersection curves still have (1,0) homology in this case, they are 
nearly touching. A reconnection does occur in Fig. IHl^d) as v is increased, causing the homology to 
become (3, 1), corresponding to a three-armed spiral. 

The intersections of the stable and unstable manifolds for when e <C 1 can be studied using 
the Melnikov technique |LM03j . The Melnikov function, M : W^{T+) M, is measure of the 
distance between the stable, V7/(71) and unstable W^{T^) manifolds to first order in e. As is 
shown in |LM03j . the Melnikov function is given by 

oo 

M= J2 dJ{Xo) o /* . (32) 

t=— oo 

Here, dJ is the one- form defined by the invariant J, H30|) . and Xq is the "perturbation vector field" 

= Pl + P2. 



Xo{x) ^ ^J.{f,\x)) 



e=0 



If M has a nondegenerate zero crossing then for e small enough there is a true transversal intersection 
of VF" and nearby. We plot contours of M in Fig. for the same set of parameters as Fig. |H1 
In this figure the domain corresponds to the fundamental domain on Wq, given by {0 > z > 
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Figure 9: Stable manifold (blue) of the circle T_ and unstable manifold (red) of the circle 7+ for the map 13111 . The 
top two panels have parameters v — 0.35 and e = 0.25 with (a) uj = 7~*, and (b) u — 7~'^. The bottom two panels 
have Lo — with (c) v = 0.30 and e — 0.35, and (d) v — 0.55, and e = 0.4..eps 
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hu{0) , < 9 < 1}, and the zero level set of M corresponds to the heavy curves. The identification 
of the upper boundary with the lower boundary to form a torus is done by shifting these boundaries 
by uj, as indicated by the arrows. Note that the homology types of the zero level sets are the same 
as those observed in Fig. [HJ 




Figure 10: Contour plots of the Melnikov function l|82^ on a fundamental domain of Wo{T+) for the same set of 
parameter values as Fig. |5| In these figures the zero level set is the heavy (black) curve, while positive and negative 
levels of M are colored red and cyan, respectively, .eps 

5 Conclusion 

The linearization of a mapping near an invariant torus with an incommensurate rotation vector 
gives rise to a quasiperiodic skew-product of the form (jlj with a periodic matrix A(6). If one can 
find a coordinate transformation in which the skew-product is rewritten as a constant matrix, then 
the system is called reducible. The eigenvalues and invariant directions of the A(9) correspond to 
the linear invariant subspaces of the torus. Because these provide the first order approximations 
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to the invariant manifolds for the invariant torus, these methods can be used in computing the 
manifolds of reducible invariant tori. We discussed the eigenvalue problem H12|l associated with the 
reducing transformation, and showed how to define generalized eigenfunctions to reduce the system 
when it has multiple eigenvalues. 

We applied our methods to compute the invariant manifolds for several examples, including 
a system studied in |LM03j by Melnikov methods. Our computations confirm the bifurcations in 
homology of the heteroclinic intersections between a pair of invariant circles predicted by |LM03j . 

While Fourier series methods have been previously used in numerical algorithms to compute 
the eigenfunctions, our simple iterative technique based on the power method works well when 
the eigenspaces are one dimensional. Intermediate eigenfunctions can also be found using iterative 
deflation methods, though their numerical stability properties — inherited from the constant matrix 
case — make them problematic. In the future we will present methods based on an iterative QR 
decomposition which much more effectively deflate A{6). Algorithms for the cases of complex and 
multiple eigenvalues will also appear in that work. 
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